%% Navier-Stokes equation solver
%
%
%

%% Problem settings

clear all
clc

addpath(genpath('./'));

global NIt
global NJt
global NIe
global NJe

%% Convergence parameters--------------------------------------------------
convTol = 10^-5;
%--------------------------------------------------------------------------

% Domain dimensions

Lxt  = 12;
Lxe  = 4;       % Longitude before expansion
Lyt  = 2;
Lye  = 1;

% Discretization, pressure finite volumes

% I denotes x direction
% J denotes y direction

elemsPerLenghtUnit = 7;

NIt = Lxt*elemsPerLenghtUnit;
NIe = Lxe*elemsPerLenghtUnit;
NJt = Lyt*elemsPerLenghtUnit;
NJe = Lye*elemsPerLenghtUnit;

% Nodal spacing

h   = Lxt/NIt;

% Constitutive parameters

nu  = 0.01;
rho = 1;

% Boundary conditions


% Time domain and discretization
dt = 0.1; % las soluciones aneriores eran con dt=0.05
T0 = 0;
T1 = 100.;
nSkip = 100;             % skips saving solution
nTimeSteps = floor((T1-T0)/dt);
nSaveSteps =  floor((T1-T0)/(dt*nSkip));
nSave      = 0;        % iniciatlization nSave 

% dtLim1 = 0.25 * h^2 / nu
% dtLim2 = 2*nu/1
% 
% if (dt>dtLim1 || dt>dtLim2 )
% fprintf ('stability warning')
% end



%% Pre-processing
%  Memory allocation

dimP = NIt*NJt-NIe*(NJt-NJe);
dimU = dimP + NJt;
dimV = dimP + NIt;

% Problem variables
Phy = zeros(dimP,1);
U = zeros(dimU,1);
V = zeros(dimV,1);

% Auxiliars: midstep variables & Pk for iterative pressure scheme computation
Pk = zeros(dimP,1);
Us = zeros(dimU,1);         % U star
Vs = zeros(dimV,1);
Ua = zeros(dimU,1);         % Previous step solution
Va = zeros(dimV,1);

% Saving solution over time
PPT = zeros(dimP,nSaveSteps);
PT = zeros(dimP,nSaveSteps);
UT = zeros(dimU,nSaveSteps);
VT = zeros(dimV,nSaveSteps);

% Auxiliar linear matrix & load vector for pressure implicit equation
M = zeros (dimP,dimP);
l = zeros (dimP,1);


% Dimensions
NIut = NIt+1;
NJut = NJt;
NIue = NIe;
NJue = NJe;

NIvt = NIt;
NJvt = NJt+1;
NIve = NIe;
NJve = NJe+1;

% Mapping spatial points

Xp = PspatialMap( h , Lxt , Lxe , Lye , Lyt , NIt , NJt , NJe , NIe , dimP);
Xu = UspatialMap( h , Lxt , Lxe , Lye , Lyt , NIut , NJut , NJue , NIue , dimU);
Xv = VspatialMap( h , Lxt , Lxe , Lye , Lyt , NIvt , NJvt , NJve , NIve , dimV);

assemblePressureMatrix;

% Vectorial implementation
%assembleDMatrix;
assembleDMatrix_zeroGradient;
Dp = sparse( - (nu/2) * D + (1/dt)* eye(dimU+dimV));
dirichletConditionsDp;
Dpp = sparse( (nu/2) * D + (1/dt)* eye(dimU+dimV)); %  <--- this is the theoretical

% A_h( u^{n})
Au_vec = zeros(dimU+dimV,1);
% A_h( u^{n-1})
Aum1_vec = zeros(dimU+dimV,1);
%% Processing - Predictor step
%
%
%
%

% Temporal loop

time = T0;
nstep = 0;
Um1 = [U; V]; % vector [u v] at previous time
aux = [U; V];

while (time<T1)
    tic
    %if (nstep>1)
    %end      
    
    time = time + dt;
    nstep = nstep+1;
    
    fprintf('Solving timestep %d of %d \n',nstep,nTimeSteps)
    
       
    % Velocity in X
    for i=1:NIut
        for j=1:NJut
            
            if (i<=NIue && j>NJue) continue; end
            
            indexU_cen = UmapOpt(i, j, NIue, NJue, NIut);
            
            U_cen   = U ( indexU_cen );  % center of finite volume
            
            if (j<=NJue)
                % Boundary conditions on inlet & outlet
                if (i==1) Us ( indexU_cen ) = getParabola ( Xp ( indexU_cen,2 ) ) ; continue ; end    % Boundary condition inlet
                if (i==NIut) Us ( indexU_cen ) = Us ( UmapOpt(i-1, j, NIue, NJue, NIut) ) ; continue ; end    % Boundary condition oulet - null normal gradient
            else
                % Boundary conditions on expansion & outlet
                if (i==NIue+1) Us ( indexU_cen ) = 0 ; continue ; end    % Boundary condition expansion    
                if (i==NIut) Us ( indexU_cen ) = Us ( UmapOpt(i-1, j, NIue, NJue, NIut) ) ; continue ; end    % Boundary condition oulet - null normal gradient
            end
            
            % vectorial ---------------------------------------------------            
%            Au_vec(Umap(i,j)) = getAx( i,j,U,V,NJue,NIue ,NJut,h);
            Au_vec(indexU_cen) = getAxOpt( i,j,U,V,NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt ,h);
            
%            Aum1_vec(Umap(i,j)) = getAx( i,j,Um1(1:dimU),Um1(dimU+1:end),NJue,NIue ,NJut,h);
            Aum1_vec(indexU_cen) = getAxOpt( i,j,Um1(1:dimU),Um1(dimU+1:end),NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt ,h);
            
        end
    end
    
    
    % Velocity in Y
    for i=1:NIvt
        for j=1:NJvt
            
            if (i<=NIve && j>NJve) continue; end
            
            indexV_cen =  VmapOpt(i, j , NIve, NJve, NIvt);
            
            V_cen  = V ( indexV_cen );  % center of finite volume
            
            if (j<=NJve)
                % Boundary conditions on top and bot boundaries
                if (j==1) Vs (indexV_cen) = 0; continue ; end
                if (j==NJve && i<=NIve ) Vs (indexV_cen) = 0; continue ; end
            else
                % Boundary conditions on top and bot boundaries
                if (j==NJvt ) Vs (indexV_cen) = 0; continue ; end
            end
            
            % vectorial ---------------------------------------------------
%            Au_vec(dimU+Vmap(i,j)) = getAy( i,j,U,V,NJve,NIve ,NIvt,h);
            Au_vec(dimU+indexV_cen) = getAyOpt( i,j,U,V,NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt ,h);
%            Aum1_vec(dimU+Vmap(i,j)) = getAy( i, j, Um1(1:dimU), Um1(dimU+1:end), NJve, NIve, NIvt, h);
            Aum1_vec(dimU+indexV_cen) = getAyOpt( i, j, Um1(1:dimU), Um1(dimU+1:end),NIue, NJue, NIut, NJut, NIve, NJve, NIvt, NJvt ,h);
        end
    end
    
    
    % Mount vector
    U_vec = [U; V];
    
    % RHS = (3/2)*(Dt) A_h( u^{n}) - (1/2)*(Dt) A_h( u^{n-1}) - D_h(u^n)
    % - (2\nu)/(Dt) u^n
    %   First step of the Adam-Bashforth got to be done by an 1 order
    %   scheme
    if(nstep ==1)
        RHS = -(3/2) * Au_vec + (1/2) * Aum1_vec + Dpp * U_vec;
    else
        RHS = -(3/2) * Au_vec + (1/2) * Aum1_vec + Dpp * U_vec;
    end
    
    % RHS = BOUNDARY CONDITIONS
    RHSboundaries;
    
    Ustar_vec = Dp \ RHS;
    
    Us = Ustar_vec(1:dimU);
    Vs = Ustar_vec(dimU+1:end);
    
    
    %% Processing - pressure step
    %
    
    c2 = h/dt;
    
    for i=1:NIt
        for j=1:NJt
            
            if (i<=NIe && j>NJe) continue; end
            
            U_right = Us ( UmapOpt(i+1, j, NIue, NJue, NIut) );
            U_left  = Us ( UmapOpt(i, j, NIue, NJue, NIut) );
            V_top   = Vs ( VmapOpt(i, j+1 , NIve, NJve, NIvt));
            V_bot   = Vs ( VmapOpt(i, j , NIve, NJve, NIvt));
            
            if (i~=NIt)
                l ( PmapOpt(i, j, NIe, NJe, NIt) ) = U_right - U_left + V_top - V_bot;
            else
                l ( PmapOpt(i, j, NIe, NJe, NIt) ) = 0 ;
            end
            
        end
    end
    
    l = l*c2;
    
    Phy = M\l;
    
    
    
    %% Processing - projection step
    
    c1 = dt/h;
    
    % Velocity X
    
    for i=1:NIut
        for j=1:NJut
            if (i<=NIue && j>NJue) continue; end
            
            indexU_cen = UmapOpt(i, j, NIue, NJue, NIut);
                        
            % Right boundary
            if ( i==NIut ) 
                U(indexU_cen) = U(UmapOpt(i-1, j, NIue, NJue, NIut)); 
                continue; 
            end
            
            % Left boundary
            if ( i==1 ) 
                U(indexU_cen) = Us(indexU_cen); 
                continue; 
            end
            if ( i==NIue+1 && j>NJue ) 
                U(indexU_cen) = Us(indexU_cen); 
                continue; 
            end
            
            U(indexU_cen) = Us(indexU_cen) - c1 * ( Phy(PmapOpt(i, j, NIe, NJe, NIt)) - Phy(PmapOpt(i-1, j, NIe, NJe, NIt)) );
        end
    end
    
    % Velocity Y
    
    for i=1:NIvt
        for j=1:NJvt
            if (i<=NIue && j>NJue) continue; end
            
            indexV_cen = VmapOpt(i, j , NIve, NJve, NIvt);
            
            % Top boundaries
            if ( j==NJvt ) 
                V(indexV_cen) = Vs(indexV_cen); 
                continue; 
            end
            if ( (j==NJve) && (i<=NIve) ) 
                V(indexV_cen) = Vs (indexV_cen); 
                continue; 
            end
            
            % Bot boundary
            if ( j==1 ) 
                V(indexV_cen) = Vs(indexV_cen); 
                continue; 
            end
                        
            V(indexV_cen) = Vs(indexV_cen) - c1 * (Phy(PmapOpt(i, j, NIe, NJe, NIt)) - Phy(PmapOpt(i, j-1, NIe, NJe, NIt)) );
        end
    end
    
    % Saving results        
    if ( mod(nstep-1,nSkip) == 0 )
        %% Calculates the pressure field
%         P = getPfromPhy( Phy, U_vec, Ustar_vec, nu, h, NIt, NJt, NIe, NJe, NIut, NJut, NIvt, NJvt, NIue, NJue, NIve, NJve, dimU, dimV, D );

        
        nSave = nSave+1;
        PT (:,nSave)= Phy;
%         PPT (:,nSave)= P;
        UT (:,nSave)= U;
        VT (:,nSave)= V;
        
    end
    
    %% convergence check---------------------------------------------------
    if (nstep > 1)
        UDiffSum = 0;
        for i = 1 : dimU+dimV;
            UDiff = abs(U_vec(i) - Um1(i));
            if (U_vec(i) ~= 0)
                UDiff = abs( (U_vec(i) - Um1(i)) ./ U_vec(i) );
            end;
            UDiffSum =  UDiffSum + UDiff;
        end; 
        err = UDiffSum/ (dimU +dimV);
        fprintf('Stationary control: %f \n', err);
        errorH(nstep) = err;
        if (err <= convTol)
            fprintf('Stationary state reached at iteration %d, convergence error: %f', nstep, err);
            convFlag = true;
            break;
        end;
    end;
    %----------------------------------------------------------------------
    
    if(nstep~=1)
        Um1 = aux; % vector [u v] at previous time
    end
    aux = U_vec; % vector [u v] at previous time
    
    toc
end


%% Post-Process

plotInTime;

